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ABSTRACT 



The standard Grad-Shafranov equation for axisymmetric toroidal plasma 
equilibrium is customary expressed in cylindrical coordinates with toroidal con- 
tours, and through which benchmark equilibria are solved. An alternative ap- 
proach to cast the Grad-Shafranov equation in spherical coordinates is presented. 
This equation, in spherical coordinates, is examined for toroidal solutions to de- 
scribe low (3 Solovev and high (3 plasma equilibria in terms of elementary func- 
tions. 



Subject headings: Toroidal Equilibria 
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I. Grad-Shafranov Equation 



The Grad-Shafranov equation (1,2) for toroidal plasma equilibria is traditionally 
formulated in a toroidal system using cylindrical coordinates to suit the plasma topology. 
This equation describes the poloidal magnetic flux, in cylindrical coordinates, with source 
terms that are themselves flux gradients. By specifying the source functionals in magnetic 
flux, toroidal equilibria are solved by different methods. Because of the cylindrical 
representation of a toroidal topology, the solutions are not easy to visualize and the 
mathematical analytic tools are rather limilted, which weaken the communication between 
theorists and experimentalists. Here, we propose to solve axisymmetric plasma equilibria 
in standard spherical coordiantes by looking for solutions with toroidal topology. With 
the familiar spherical system, and the large volume of well-known special functions in this 
coordinate system, the solutions in toroidal topology appear in a much friendly presentation, 
which improves physical interpretations. 

We first derive the Grad-Shafranov equation in spherical coordinates by considering 



With axisymmetry, the magnetic field and the current density can be represented by two 
scalar functions in standard spherical coordinates 



J X S- Vp = , 
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Making use of Eq. [2] to eliminate the current density, Eq. [T] renders three components. The 
(p component contains only the magnetic force, and it is 



dr de de dr ~ ' ^ ' 

Qir,e) = Q{P{r,e)) . (5b) 



As for the 9 component, it reads 

Since the spatial dependencies of Q{r, 0) is through a functional of P(r, 6), we could express 
the 9 dependency in plasma pressure by P to write p(r, 9) = p{r, P), and the above equation 
becomes 



d'P 1 d 1 dP dQ 2 . 2ndp 

{^ + ;:^^-^^(-^^) + ^^} + /^- 

This equation is the Grad-Shafranov equation (1,2) of axisymmetric toroidal plasma 
equilibrium, represented in spherical coordinate system. The first three terms represent the 
nonlinear force-free field with fiJ = K{P)B, where K{P) = dQ/dP is a scalar function. 
This can be verified from Eq. [3] and Eq. Hlwhen we impose fiJ = K{P)B. In particular, we 
would have the linear force- free field should we take Q{P) = aP with constant K{P) = a. 
The last term is the plasma pressure balance. The magnetic function Q'^{P) and the 
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pressure function p{P) are source functions of the Grad-Shafranov equation that need to be 
specified. Finally, the r component of Eq. [T] reads 



p]- ^ .29Pfd^p 1 d 1 dP idQ\ 

dr^^^' ^~ /iVsin^^ dr^dr^ r^^^^ dO^smO d9 ^ ^ 2 dP ^ 

dP d 

The term dp{r, P)/dr on the left side refers to the total radial derivative, which has an 
explicit and an implicit part, and the right side of this radial equation is the magnetic force. 
Making use of the 9 component, Eq. [3, the right side is equal to the implicit part of the 
radial pressure gradient, which cancels the same term on the left side, the implicit part of 
the radial derivative. The radial component, therefore, reads dp/dr = 0, and the plasma 
pressure is a function of P only. 



p = p{P) 



(9) 
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II. Low (3 Equilibria 

To solve Eq. [7] analytically, we take the source functions as 



Q\P) = 2b^P + a^P^ + Ql 
p{P) = poP + p'^lnP-Co 



(10a) 
(10b) 



where po and p'^ carry the dimension of pressure, and 6^ measure the dimension of 
magnetic field, Ql and Cq are constants, and P(r, 9) is a dimensionless function. With 

= and pg = 0, Solovev (3) solved the Grad-Shafranov equation for low /? equilibria with 
linear source functions. These source functions have been extended by other investigators 
(4,5). With ^ 0, the quadratic form for plus a hkewise quadratic form for p, without 
the InP term, have been analyzed by many authors in different ways (6-9). In particular, 
Hu has considered a cubic form for (10). On the other hand, numerical routines have 
also been advanced by many authors (11,12). Here, we include a InP term in the plasma 
pressure function, and regard Eq. llOal and Eq. llObI as the generalized Solovev case, and 
solve for the Grad-Shafranov equation in spherical coordinates for a low /3 plasma. The Q 
component, therefore, becomes 

d^P d 1 dP 1 

{r^— - + sin ^— (—-—) + (ar)2p+(6r)2} + /iporSin2^ + wSin2^- = . (11) 
ar^ od smfe* oO P 

Dividing this equation over by P, separating the variables by writing P{r,6) = R{r)<d{6), 
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denoting x = cos6, and with n{n + 1) as the separation constant, we have 



r'-^^ - n{„ + l)R{r) = -(ar)=fl(r) - ^|5(6r)-'<^^^ . (12b) 



These equations can be solved by iterations, because of the low /5 plasma. We, therefore, 
first solve the homogeneous equations with the right sides null giving 



Qoix) = (1-^')^^ , (13a) 
Ro{r) = Ao(6r)"+i , (13b) 



where Pn{x) is the Legendre polynomial. The coefficient Ao reflects the amplitude of the 
magnetic flux through Eq. [3l To iterate on the zeroth order solutions, we take n = 1 and 
substitute Eq. I13al and Eq. llSbl to the right sides of Eq. I12al and Eq. Il2bl to get 



{brY ^p'q (br)^ (1 — x^) 



n{n + l)Ri{r) 



(ar?Rn(r) 
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(14b) 
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By simple inspection, the first order solutions and the general solutions are 



ei(x) -- 

Ri{r) = - 



2A* 



10^ ^ 



R{r) = Aoibrf - ^(6r)^ 



(15a) 



(15b) 



We see that G(x) has a lobed solution centered at x = with a maximum amplitude 
(1 — 1/2Aq). The radial solution R{r) is composed of two terms. The first term is positive 
and increases to the second power of (6r), while the second term is negative and increases 
to the fourth power of (br). The radial function R{r) is positive at small (br) and negative 
at large (br). It vanishes at (6r)^ = and (6r)^ = IOAq/Ai, with a maximum in between. 
This maximum of R{r) and the lobed solution of 0(x) are essential for a torus structure for 
the magnetic flux. With the given source functions, we notice that the iteration scheme 
can not be carried beyond the first round. Consequently, if Eq. I15al and Eq. I15bl have any 
meaning at all, Ri{r) has to be much less than Ro{r), and 0i(x) has to be much less than 
6o(x), leading to Ai/Aq « 1 and 1/Aq « 1. This requires [(a/6)^ + (fipo/b'^Ao)] « 1, 
or {a/by « 1 and {^po/b^Ao) « 1, and [1 + (fipQ/b^Ao)] « Aq, for a low /3 plasma. 
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III. Magnetic Torus 

With the spatial structure solved, the magnetic field components are given by 



1 IdP 
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(16a) 
(16b) 
(16c) 



The solution R{r) vanishes at some r where we have Br{r) — 0. The solution Q{x) also 
vanishes at some x. Together they describe the magnetic fields. Within this region of (r, x), 
the topological center defined by dR{r)/dr = and dQ{x)/dx = has Br = and Bq = 0. 
This is the magnetic axis, r — r^, where the magnetic field is entirely toroidal. The field 
lines about this center are given by 



Br _ Be _ 
dr rd9 

By axisymmetry, the third group is decoupled 
lines on an (r — 9) plane, we consider the first 
P — R{r)Q{x) equals to a constant or 



r sin 9d(f) 

from the first two groups. For the field 
equality between Br and Be which gives 



P{r,x) = R{r)e{x) = C . (18) 

In other words, the nested poloidal field lines are given by the contours of P{r,x) on the 
(r — x) plane. At the topological center, we have Q{x) maximum and R{r) maximum, so 
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that P(r, x) is maximum. Since r sin 6 is the distance of a point on the (r — x) plane to the 
z axis, Eq. I16cl states that the hne integral of around the circle on the azimuthal plane 
is measured by 27!-Q, 

27rr sin 6B^ = 2'kQ = ^1^ , 

and the center has the maximum of this line integral about the axis of symmetry. Also, it is 
evident that Q is equivalent to the axial current, where the constant part Q = Qo amounts 
to a uniform current. As for P, we evaluate the poloidal magnetic flux by integrating 
Eq. I16bl on the x = plane over a cross section to give 




2nrBedr = -2n{P{z) - P{z,)) . 
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IV. Near Nonlinear Force-Free High [3 Equilibria 

For high (3 plasmas, we take source functions 



Q\P) = a^P^ + Ql , (19a) 
p{P) = poP - Co . (19b) 



The Grad-Shafranov equation, Eq. [3, reads 



y^ + sme-^i--^) + iaryP} + f^Por'sm'e = . (20) 

The first three terms correspond to the force-free magnetic field as can be verified by setting 
/X J = aB between Eq. |3] and Eq. HI and the last term is the plasma pressure balance term. 
Separating the variables gives 



X 



+ [(ar)^ -n{n + l)]R 



+ n{n+ l)0(x) = 

1^ 



a" 



e 



(21a) 
(21b) 



The first equation gives 



e(x) = (1-^^)^^ = (1-^^) , (22) 
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where n = 1 is again chosen. This solution gives 0(x) positive for any x, which is 
appropriate for small aspect ratio plasmas. As for the second equation, having n = 1, the 
solution is given by R{r) = Ro{r) +Ri{r), where -Ro(^) and Ri{r) are the homogeneous and 
particular solutions. The homogeneous solution is described by 



Ro{r) = Aoarjn{ar) = Aozjn{z) , (23) 

where jn{z) is the spherical Bessel function. As for the particular solution, denoting 

Ai = {fipo/a'^), we have 



z^^^ + [z^ -n{n + l)]Ri = -A^z^ , 

i?i(r) = -Ai{arf = -A^z^ . (24) 

We note that the homogeneous solutions, Ro{r) and Q{x), correspond to the linear or 
nonlinear force-free solutions of the first three terms of Eq. [20] with null or finite Qq 
respectively. The last term for plasma pressure appears only in the particular solution, 
i?i(r), that keeps the pressure balance. The homogeneous radial solution is an oscillating 
function in z = ar, which has sucessive maxima, and the homogeneous meridian solution 
has a lobe peaked at x = 0, like the low (3 case. The superposition of the particular radial 
solution only slightly modifies the homogeneous solutions. We could use the region between 
z = and the first root of jn{z), with n = 1, to describe low aspect ratio high (3 toroidal 
plasma equilibria. Since Eq. [22], Eq. [23l and Eq. [21] are exact solutions of Eq. [201 with 
source functions of Eq. I19al and Eq. Il9bl these exact solutions decribe high (3 plasmas, that 
are nearly force-free. It is interesting to note that exact near force-free solutions appear in 
high P equilibria, while the low (3 equilibria are not force-free. 
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V. Low Aspect Ratio Fusion Tori 

These benchmark analytic equihbria are most relevant for current international 
collaboration scale fusion tokamaks, such as JET and ITER, with D shaped plasma 
cross section. The homogeneous radial force-free solution Ro{r) is shown in Fig.l with 
Aq = 0.05 Tm^ taken arbitrarily. This Ro{r) solution has the first null point at zq = 4.5. 
To consider the particular solution Ri{r), we first have to determine the constant 
Ai = {fipo/a^). To fix the parameter a, we remark that the first stationary point, for the 
magnetic axis, is around 

Taking the major radius = 0.5 m for a given low aspect ratio spherical tokamak, 
we would have a = z^jr^ = 5A/m. Considering plasma pressure of 10^ Pa, we have 
Ai = 1.3 X 10~Vo = 1-3 X 10~^. With the given Ai, the particular solution of Eq. [2H is 
shown in Fig.2. Superimposing the two solutions gives the near force-free general solution 
R{z) = Ro{z) + Ri{z) which shifts the null point of R{z) to zq < 4.5, where the plasma 
pressure vanishes. The corresponding contour point of Zq < 4.5 on the inner side of the 
torus is 2;o = 0, where the plasma pressure also vanishes. However, there are engineering 
aspects of fusion torus, such as divertor scrape-off effect, that alter the outer and inner zq 
plasma boundaries. This effectively means that the plasma boundary is set by the constant 
Co in Eq. Il9bl or Eq. llObI determined by engineering considerations. 

The poloidal magnetic contours of Eq. [T8]are shown in Fig. 3. Since the plasma pressure 
of Eq. llQbl is expressed in terms of P, it shares the same D shaped contours of Fig. 3. As for 
the toroidal field, with Qq = 0, the toroidal field circulation about the z axis is expressed in 
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terms of 2ttQ = 2TTaP, and therefore it also shares the contours of Fig. 3. As for the current 
density of Eq. HI making use of the Grad-Shafranov equation of Eq. 120] gives 



rsmtJ r do or 



{+i^,_^,Wis,e} . (25) 



rsin6' r 89 dr 

The middle equality is obtained by cancelling the Ri term with the ^po term using Eq. [2 
Analogous to the magnetic field lines, the current density field lines are given by 



7^ = 4 = — ^ ■ (26) 
dr rdU rsmf 



Considering the first equality, the poloidal current density contours are given by 



Q{r,x) = C . (27) 

As a result of this. Fig. 3 also represents the poloidal current density contours in the linear 
force-free case of Ql = such that Q{P) = aP{r, x) = aR{r)Q{x). 

The field along x = cos ^ = in Tesla, Eq. I16cl is shown in Fig. 4. Since this field 
scales as R{r)/r, it has a ji(r) profile of Fig.l, and it peaks around z = 2.1. The Bg field 
along X = cos 6 = in Tesla, Eq. I16bt is also shown in Fig.4. This component scales as 
{l/r){dR{r)/dr), and it crosses zero around z = z^, = 2.6 which sets the magnetic axis. We 
note that the magnetic axis at z = 2.6, where Bg vanishes, and the stationary point of 
the toroidal field at z = 2.1, where B^f, is maximum, do not coincide. They are separated 
by Az = aAr = 0.5, or Ar = 0.1m. Since the toroidal field contours are off-centered 
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from the poloidal field contours, the magnetic pressure contours are not aligned with the 
plasma pressure contours. The discrepancy is due to the magnetic curvature which gives 
rise to magnetic tension according io jiJ x B = [B ■ V)B — V-B^/2. This difference A2; 
would be reduced if we use subsequent second and third peaks of RQ{r) — 6.1,9.3 and 
corresponding ji(r) = 6.0,9.2, not shown in Fig.l. As an example, for the second peak, 
we would have = 6.1 and A2; — 0.1. With major radius — 2m, this would give 
Ar = {/S.z/z^)r^ — 0.03 m, and would have high aspect ratio equilibria. The toroidal and 
poloidal current densities, and Jq, along x = cos 9 = in 10^ A/ are shown in Fig. 5. 
Besides the scale difference, Fig.4 and Fig. 5 are slightly different with J^/B^f, slightly larger 
than Jb/Bq. Had not been the plasma pressure dependent particular solution Ri{z), they 

— * — * 

would be the same, because of the force-free nature with J and B being proportional. 
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VI. Safety Factor 



One important parameter that measures the quahty of equihbrium is the safety factor, 
which measures the number of turns that the toroidal field circles about the toroidal axis of 
symmetry through the major radius as the poloidal field makes one turn about the magnetic 
axis through the minor radius returning to the same poloidal angle, but not necessarily the 
same toroidal angle. In spherical coordinates, this amounts to the winding number of the 
toroidal field about the z axis of symmetry as the poloidal field circles one complete turn 
about the magnetic axis. With our spherical coordinates, the poloidal field is composed of 
Br and Bg, and the usual safety factor expression does not apply. In order to rederive an 
expression for the safety factor, we consider the field line equation, Eq. [171 Since Br and 
Bg are constrained by the poloidal field line contours, we need only one of them to track B^ 
around the torus. By considering the radial component, and writing Qq = a^7^, such that 
Q{P) = a{P'^ + 7^)^''^, the field line around the torus is given by 



where the integration is carried out over the contour R{z)Q{x) = C, as in Fig. 3, with x{z) 
in the integrant given by 

The integration is bounded hj Zi < z < Z2 on the a; = horizontal plane, where 




(28) 
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Xi = x{zi) = and X2 = x{z2) = are the first and second intercepts of a given contour 
on tlie liorizontal axis. Altliougli tlie integrant is singular at tlie bounds zi and Z2, it is 
integrable. As tlie contour value C increases towards the maximum of R{z), the contour 
closes onto the magnetic axis of z = z^,, and xi and X2 come closer together towards z*. In 
this limit, x{z) in the denominator approaches zero, and the integration interval zi and 
Z2 of Eq. [28] also reduces to null. As a result, the contour integral near the magnetic axis 
depends on the limiting integral of the last factor dz/x{z). To get the correct limit, we 
Taylor expand R{z) about z^ to second order to obtain 



[R\z.)+j']y' P dz_ _ , 2[R\z.)+^^] \z2-zA 
" R{z.) I, x{z)~ ^R{z.)\dm{z.)ldz^\^ \z,-zA- ^""^^ 

Since the R{z) profile is not symmetric around 2;^,, as shown in Fig.l, zi and Z2 of a given 
contour value do not approach z^ at equal intervals, and 2TTq{zt,) is finite. Since the second 
order term of the Taylor expansion is a symmetric quadratic term in [z — z*), this left-right 
non-symmetry of the maximum comes from the third order cubic term. 

As an example, let us consider 7^ = 0.01. Since the maximum value of P = RQ is 
about 0.05, according to Fig.l with G = 1, the constant 7^ term dominates the term 
in Q. By Eq. I19al this amounts to an almost uniform axial current /i/^ = 2'kQ. In this 
case, the scalar K{P) = dQ/dP is not a constant, but a function of P, which makes the 
force-free equation nonlinear. The difference between the linear and nonlinear force-free 
fields is that the toroidal magnetic field vanishes at some radial positions for linear case, 
while it is always poisitive definite for nonlinear case. With the corresponding parameters 
in Figs. 1-5, the safety factor profile is shown in Fig. 6, which ressembles well with laboratory 
profiles, having a basin within the plasma interior surrounded by high values on the border. 
The value at the border can be lowered by reducing the constant 7^. 
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VII. Conclusions 

We have examined the axisymmetric toroidal plasma equilibria for low and high (5 
plasmas. The center of the spherical system, z = 0, has been excluded because this point 
usually lies outside the plasma domain due to the machine vessel. However, for spheromaks, 
the center is connected to the plasma domain. In this case, our solutions provide hollow 
profiles for plasma pressure. Rather than solving the Grad-Shafranov equation in toroidal 
geometry with cylindrical coordinates, devised to envelope axisymmetric toroidal plasmas, 
we have taken a different approach to represent the Grad-Shafranov equation in spherical 
coordinates and have solved for toroidal solutions. This approach allows benchmark plasma 
equilibria, such as the generalized Solovev type low /3 equilibrium and the standard high (3 
equilibrium, be solved by separation of variables in terms of elementary special functions, 
and provides easier visuahzation of the solutions. The essential features of the solutions in 
spherical coordinates are that the radial function R{z) has a maximum aX, z — ^ and 
the meridian function Q{x) has a lobe at x = 0. These two features assure the function 
P{z,x) = R{z)Q{x) an axisymmetric structure of a torus. Reahstic high (3 near force- free 
D shaped contours are evident in these analytic representations. 

In the generalized Solovev type source functions, the equilibrium is solved approximately 
by sucessive approximation to the first order. As a result, the validity of the solution 
requires strictly a low (5 plasma. For more specific source functions, the equilibrium is solved 
exactly in closed form. The solution is, therefore, valid for high j3 plasmas. Surprisingly, the 
high (3 exact solution is constructed on a force-free magnetic field superimposed by a high 
(5 plasma contribution. Although with high /3, the plasma contribution is small comparing 
to the force-free part, making the overall solution near force-free. The reason for this near 
force-free solution to confine plasmas is that it relies basically on the tensional force of 
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the magnetic field line curvature, not the scalar magnetic pressure, to balance the plasma 
pressure. This curvature is evident because the magnetic axis where Bg = is significantly 
displaced from the maximum of B^. For this reason, the magnetic field profile remains near 
force-free, while plasmas are confined. Such strong curvature is inherent from low aspect 
ratio plasmas. The absence of such curvature in high aspect ratio plasmas obliges the 
machines to operate in low /? limit with magnetic fields away from force-free configuration. 
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Fig. 1. — The function Rq{z) = A^zji^z) with Aq = 
solution is plotted function of z in thick line, 
also plotted in thin line for comparisons. 



0.05 of the force-free homogeneous radial 
The spherical Bessel function AqJi^z) is 




Fig. 2. — The function Ri{z) = —Az^ with po = 10^ Pa and A = 1.3 x 10 ^ of the particular 
radial solution is plotted as a function of z to show the plasma pressure effect. 
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4- 




z(x=0) 



Fig. 3. — The magnetic field lines for spherical tokamak are shown with increasing contour 
labels from outer to inner contours. These contours are shared by the plasma pressrue and 
the toroidal field line integral. The axes are labelled in z together with x = cos 6. 




Fig. 4. — The toroidal and the poloidal magnetic fields, i?,^ and Bg, with Q = aP and 
are shown as a function of z along x = 0. 
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Fig. 5. — The toroidal and the poloidal current densities, and Je, with Q = aP and 
a = 5.4 are shown as a function of z along x = in units of 10^ A/m^. 
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Fig. 6. — The safety factor profile of the plasma equilibrium with nonlinear force-free homo- 
geneous solution, 7^ = 0.01, is shown as a function of z along a; = 0. 



